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Site diluted spin-1/2 Ising and spin-1 Blume Capel (BC) models in the presence of transverse field 
interactions are examined by introducing an effective-field approximation that takes into account 
the multi-site correlations in the cluster of a considered lattice with an improved configurational 
averaging technique. The critical concentration below which the transition temperature reduces to 
zero is determined for both models, and the estimated values are compared with those obtained by 
the other methods in the literature. It is found that diluting the lattice sites by non magnetic atoms 
may cause some drastic changes on some of the characteristic features of the model. Particular 
attention has been paid on the global phase diagrams of a spin-1 BC model, and it has also been 
shown that the conditions for the occurrence of a second order reentrance in the system is rather 
complicated, since the existence or extinction of reentrance is rather sensitive to the competing 
effects between D/J, fl/ J and c. 
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I. INTRODUCTION 

Investigation of disorder effects on the critical phenomena has a long history and there have been a great many of 
theoretical studies focused on disordered magnetic materials with quenched randomness where the random variables 
of a magnetic system such as random fields [j], [|J or random bonds [!, 0] may not change its value over time. On the 
other hand, site diluted ferromagnets constitute another example of magnetic systems with quenched disorder such 
as a compound A x Bi_ x C where magnetic A atoms in a pure magnet AC are replaced by non-magnetic B impurities. 
Formerly, Sato et al. [5] have shown that in a dilute lattice a Curie or a Neel temperature does not appear until 
a finite concentration of magnetic atoms is obtained if the atomic distribution is random. They have also found 
that this concentration depends on the coordination number of the lattice. After this seminal work of Sato et al. 

much attention has been paid to site dilution problem and the situation has been handled by a wide variety of 
techniques such as Bethe-Peierls- Weiss (BPW) method @, renormalization group (RG) techniq ue B -tlfjj. correlated 
effective field theory (CEFT) (TT| - [l3| . effective field theory (EFT) based on decoupling (or Zernike |l4j ) approximation 
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(DA) fl5T - [33l ]. an integral representation method [34|, Monte Carlo (MC) simulation technique [35}-|38j], Bogoliubov 
inequality approach [39j, Bethe-Peierls approximation (BPA) [40(, finite cluster approximation (FCA) which gives 
results identical to those obtained by EFT for a one spin cluster [4lT - l44j ]. third order Matsudaira approximation 
(45j . EFT with probability distribution technique [461453 ] and cluster variational method (CVM) (53|. Among the 
theoretical works mentioned above, some of the authors extended the standard dilution problem to more complicated 
versions by taking into account the transverse field interactions random fields and random bonds, as well as 
bilinear and biquadratic exchange couplings and crystal field interactions for the systems with S > 1/2. 

Mean field theory (MFT) of site dilution problem predicts that the system always has a finite critical temperature 
and stays in a ferromagnetic state at lower temperatures, except that c — where c denotes the magnetic atom 
concentration. Therefore, it is not capable of locating a critical site concentration at which the transition temperature 
reduces to zero. The reason is due to the fact that MFT neglects single-site and multi-spin correlations. On the 
other hand, EFT based on DA accounts all the single site correlations, but it also neglects multi-spin correlations 
between different sites. Hence, EFT provides results that are superior to those obtained within the traditional MFT. 
Furthermore, CEFT which is an extension of EFT partially takes into account the effects of multi-spin correlations 
and improves the results of conventional EFT in many cases. Based on the physical aspects of the problem, whether in 
EFT or CEFT formalism, evaluation of configurational averages emerging in definition of spin identities plays a critical 
role. However, as mentioned by Tucker [231 ] . the conventional configurational averaging technique applied in Refs. 
[ill [l2l . [ill [22l . l27l [3ll [32! ] is based on a procedure that decouples the site occupation variable from the thermal average 
of spin variables, even when both quantities referred to the same site while in Refs. [lSj l23l426l 123 . 1 441 ], the authors 
used an improved configurational averaging method in which only the correlations between quantities pertaining to 
different sites are neglected (decoupled). However, it is possible to improve the accuracy of these methods by including 
multi-site, as well as single site correlations. 

In this paper, we describe a new type of EFT method for investigating the thermal and magnetic properties of a 
site diluted Ising model on 2D lattices. Recently, we have successfully applied our method to random bond [5|| and 
random field [56j problems on 2D and 3D lattices with various coordination numbers. As we emphasized in these 
previous works, an advantage of the approximation method proposed by this method is that no decoupling procedure 
is used for the higher-order correlation functions. Therefore, it is expected that the accuracy of the results obtained 
within the present work may improve those of the works based on conventional and improved DA. For this purpose, 
we organized the paper as follows: In Sec. [H] we briefly present the formulations. The results and discussions are 
presented in Sec. Mil and finally Sec. [IV] contains our conclusions. 

II. FORMULATION 

In this section, we give the formulation of the present study for site diluted spin- 1/2 Ising and spin-1 Blume Capel 
(BC) models on 2D lattices. As our model, we consider N identical spins arranged on a 2D regular lattice. Then 
we define a cluster on the lattice which consists of a central spin labeled So and q perimeter spins being the nearest 
neighbors of the central spin. The cluster consists of (q + 1) spins being independent from the value of S. The 
nearest-neighbor spins are in an effective field produced by the outer spins, which can be determined by the condition 
that the thermal average of the central spin is equal to that of its nearest-neighbor spins. In the following subsections, 
we give a detailed discussion of how the present method can be formulated for spin-1/2 Ising and spin-1 BC models 
with quenched site dilution. 

A. Site diluted spin-i system 

As a site diluted spin-1/2 Ising model, we consider the following Hamiltonian 

h = -j acjSfS?, (l) 

<i,j> 

where the summation is over the nearest-neighbor pairs of spins and the operator Sf takes the values Sf ~ ±1. We 
assume that the lattice sites are randomly diluted and Cj denotes a site occupation variable which equals to 1 if the 
site is occupied by a magnetic atom or to if it is empty. 

According to the Callen identity [57( for the spin-1/2 Ising system, the thermal average of the identity aSf at the 
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site i is given by 



c i ({/ i Kf> = c i ({/ i }tanh 



f3d jJ2 c i S i 



(2) 



where j3 — j expresses the nearest-neighbor sites of the central spin and {fi} can be any function of the Ising 

variables as long as it is not a function of the site. Applying the differential operator technique [H, [59| in Eq. @ 
and using the relation 



exp(aci) = Ci exp(a) + 1 — c,, 



with the fact that c" = we get 



Ci ({fi}Sf) = Cl ({fi} j[ exp {J Cj S?V)j tanh(/3z)| a=0 . 
By putting Eq. into Eq. we obtain 

Ci ({fi}Sf) = d /{fi} J] { Cj cosh(JV) + Ci 5? sinh(JV) + 1 - Cj }\ tanh(/3x)U =0 , 



(3) 



(4) 



(5) 



where V is a differential operator, q is the coordination number of the lattice, and (...) represents the thermal average. 
Eq. (O is valid only for a given specific magnetic atom configuration. Hence, if we consider configurational averages 
then we may rewrite Eq. (O as 

(ci({fi}S!)) r = ^^{/Onic.cosh^ + Cj^smh^ + l-c^yy tanh(#B)| s=0 , (6) 

where (...) r represents random configurational averages. When the right-hand side of Eq. ^ is expanded, the multi- 
site correlation functions appear. The simplest approximation, and one of the most frequently adopted is to decouple 
these correlations which is called decoupling approximation (DA). In conventional manner, eliminating the term Cj 
from both sides of Eq. ([5]) then performing the configurational average with {fi} — 1 leads to the following equation, 



«S?)) r = cosh ( JV ) + <hSj sinh(JV) + 1 - Cj }jj tanh(/3x)U =0 - 



(7) 



In conventional DA one expands the right-hand side of Eq. ([7]) then decouples the multi-site correlations according 
to 



with 



((d...c j c k S* k ciSi...c m S* rn )) r S (c t ) r ... {c ) r (ck) p ({S%)) r ( Cl ) r ((S[)) r ... (c m ) r ((S z m )) r 



(c a ) r = c and {{S z a )) r = m a = i, ...j,k,l, ...,m. 



(8) 



However, this approximation decouples the site occupation variable from the thermal and configurational averages of 
spin variable, even when both quantities referred to the same site. 

On the other hand, an improved version of decoupling approximation deals with the quantity (c, {Sf )) . In other 
words, in an improved decoupling procedure, one expands the right-hand side of Eq. ^ instead of Eq. ([7]) and 
decouples the multi-site correlations according to 



with 



((c l ...c J c k S z k c l S[...c m S: n )) r - (ci) r ...( Cj ) r (c k (S z k )) r (c, (Sf)) r ... (c m (5* » r 
(ct) r = (cj) r = c and (c a (S*)) r = m a = k,l,...,m 



(9) 
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In this approximation, only the correlations between quantities pertaining to different sites are neglected. A detailed 
discussion about these configurational averaging techniques is also given by Tucker |23| . Whether conventional method 
or improved one, the papers which utilize these approximations claim that if we try to treat exactly all the spin-spin 
correlations emerging on the right-hand side of Eqs. © and 0, the problem becomes mathematically intractable. In 
order to overcome this point, recently we proposed an approximation that takes into account the correlations between 
different sites in the cluster of a considered lattice [111 [56|. Namely, an advantage of the approximation method 
proposed by those studies is that no decoupling procedure is used for the higher-order correlation functions. 

We state that hereafter, we will carry on the formulation of the dilute spin-1/2 system for a honeycomb lattice 
(q = 3), however a brief explanation of the method for a square lattice (q = 4) can be found in [3] Now, if we expand 
the right-hand side of Eq. ([5]) for q = 3 without using DA, we get some certain identities in the form 

((c i ...c j c k S* k c l Sf...c m S* m )} r = (c i ...c j ) r {{c k S z k c l Sf...c m S* m )) r . (10) 

In Eq. (|10p . we use the fact that occupation number Cj of a given site i is independent from the thermal average, as 
long as the correlation function does not contain a spin variable Sf, and the site occupation numbers pertaining to 
different sites are assumed to be statistically independent from each other. Hence, we may rearrange Eq. (ITU1) as 

{(c i ...c j c k SfaSf...c m SZ l )) r = (*) r ... < C ,) r ((c k Sic l St...c m S* m )) T ■ (11) 

where (ci) r = (cj) = c. In the present formulation, it is clear that Eq. ([TT|) improves EFT based on Eqs. © and © 
by taking into account the multi-site correlations. With the help of Eq. (TTT1) . and by expanding the right-hand side 
of Eq. ^ for the central site cqSq with {/,} = 1 we have 

m = ((c So)) r =xi = (3c - 6c 2 + 3c 3 )x 4 Kx + (6c 2 - 6c 3 )x 4 K 2 + 3c 3 x 4 K 3 + cx 6 K 4 . (12) 

where the terms Xi in Eq. (|12p are defined in [A] In obtaining Eq. (TT21 we use the fact that tanh(/3x) is an odd 
function. Hence, only the odd coefficients give non-zero contribution which can be given as follows: 

Ki — sinh(JV) tanh(/3x)| K= o, 

K 2 = cosh(JV)sinh(JV)tanh(/3x)| x=0 , 

K 3 = cosh 2 (JV)smh(JV)tanh(/3a;)| x= o, 

K 4 = sinh 3 (JV)tanh(/3a;)| a:= o- (13) 

For comparison, if we apply the improved decoupling approximation given in Eq. Q then Eq. (|12[) reduces to 

m = (3c - 6c 2 + 3c 3 )mi^i + (6c 2 - 6c 3 )mK 2 + 3c 3 mK 3 + cm 3 K 4 , (14) 

which is identical to those obtained in Refs. [ltl [2^, [25[ . Additionally, applying the conventional method ([5]) gives 
the following result 

m = (3c - 6c 2 + 3c 3 )mi^i + (6c 2 - 6c 3 )mK 2 + 3c 3 mK 3 + c 3 m 3 K 4 . (15) 

It seems like it is fortuitous that although, the equations of states of approximations ((5J) and @ are differ from each 
other in the last term, they give the same phase diagram in (fcsT c /J — c) plane. The reason comes from the fact that 
both approximations ignore the term m 3 in the limit T — > T c . Hence, it should be emphasized that the importance 
and distinction of our method becomes evident by expansion of Eq. ^ without using any kind of DA. 

The next step is to carry out the configurational and thermal averages of the perimeter site in the system, and it 
is found as 

(({fs}c S S 5 )) r = (c 5 ({f s } (cq cosh( JV) + c S sinh( JV) + 1 - c )» r tanh(/3(x + 7 )). (16) 

From Eq. (fT6|) with S — {fs} = 1 we get the following identity 

«ciSi» r = x± = (c - c 2 )A 1 + c 2 A 2 + c Xl A 3 . (17) 

For the sake of simplicity, the superscript z is omitted from the left- and right-hand sides of Eqs. (fT2)) and (fT7|). The 
coefficients in Eq. (| 1 T[) are given as 

Ax = tanh(/3(.T + 7 ))| a . =Q , 

A 2 = cosh(JV)tanh(/3(a- + 7))| x=0 , 

A 3 = sinh(JV)tanh(/3(x + 7))U=o- (18) 



5 

The coefficients in Eqs. (fl~3|) and ((18)) can easily be calculated by applying a mathematical relation, e aV f(x) = f(x+a). 
In Eq. (fT8"|) . 7 = (q — l)A is the effective field produced by the (q— 1) spins outside of the cluster, and A is an unknown 
parameter to be determined self-consistently. 

Eqs. (fT2|) and (fTTj) are the fundamental correlation functions of the system. On the other hand, for a honeycomb 
lattice, taking Eqs. ([6]) and (|16p as basis, we derive a set of linear equations of the site correlation functions in the 
system. At this point, we assume that (i) the correlations depend only on the distance between the spins and (ii) the 
average values of a central site and its nearest-neighbor site (it is labeled as the perimeter site) are equal to each other 
with the fact that, in the matrix representations of spin operator S, the spin-1/2 system has the property (S) 2 = 1. 
Thus, the number of linear equations obtained for q = 3 and q = 4 reduces to six and eight, respectively, and the 
complete sets are given in [X] 

Finally, e.g. if Eq. (|A2j) for q — 3 is written in the form of 6 x 6 matrix and solved in terms of the variables Xi 
(i = 1,2,. ..6) of the linear equations, all of the site correlation functions can be easily determined as functions of 
the temperature and Hamiltonian parameters. Since the thermal and configurational average of the central site is 
equal to that of its nearest-neighbor sites within the present method, the unknown parameter A can be numerically 
determined by the relation 

x\ = X4. (19) 

By solving Eq. ((19)) numerically for a given fixed set of Hamiltonian parameters, we obtain the parameter A. Then 
we use the numerical values of A to obtain the site correlation functions which can be found from Eq. (|A2[) . Note 
that A = is always a root of Eq. (fT9|) corresponding to the disordered state of the system whereas the nonzero root 
of A in Eq. (|19p corresponds to the long-range-ordered state of the system. Once the site correlation functions have 
been evaluated then we can give the numerical results for the thermal and magnetic properties of the system. Since 
the effective field 7 is very small in the vicinity of fcsT c /J, we can obtain the critical temperature for the fixed set 
of Hamiltonian parameters by solving Eq. (fT^|) in the limit of 7 — > 0, and then we can construct the whole phase 
diagrams of the system. 



B. Site diluted spin-1 Blume-Capel model with transverse field interactions 



Site diluted spin-1 Blume-Capel (BC) 60i, 61] model with transverse field interaction is represented by the following 
Hamiltonian 

H = -JJ2 ^iSfSf -Dj2c t (S[f - £ c t Sf, (20) 

<i,j> i i 

where Sf and Sf denote the z and x components of the spin operator, respectively. The first summation in Eq. (|20[) 
is over the nearest- neighbor pairs of spins and the operator Sf takes the values Sf = 0, ±1. J, D and f2 terms stand 
for the exchange interaction, single-ion anisotropy (i.e. crystal field) and transverse field, respectively. 
By using the approximated spin correlation identities [62[ 

Ci {{ fi}Si } = (k ( {/;}— , . ), (21 

\ Tr l exp(-/3i? i ) / 



c i ({/ i }(«9?) 2 > = c^{/ i }- 



Tr t (5f) 2 cxp(-/3g. t 
Tr, exp {-(3 Hi) 



(22) 



and following the same methodology of Sec. Ill A[ we can obtain the general form of the site correlation functions for 
the central site as follows 



(ci ({fi}S*)) r 



n h ( s rf cosh( jv) + cjS i sinh( jv) + 1 ~ Cj (^l ) ) f ^ x =°-- 



(23) 



( Ci ({fi}(S?) 2 )) r = /« /{ft J] [cj (S|) 2 cosh(JV) + CjS' s sinh(JV) 



{si 



G(x)\ x=0 , 



(24) 
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where the functions F(x) ad G(x) can be found in Ref. [Hj]. By expanding the right hand side of Eqs. (l23|) and (|24|) 
according to Eq. (flT]) for c Sq and Co(S'o) 2 , respectively with {/j} = 1 and taking only the nonzero terms we get 
magnetization and quadrupolar moment of the central site as follows 

m = ((coSo)) r = xi = 3cx 4 k\ + cxek3 + (— 6A»i + 6k 2 )cxs + (3fci — 6k 2 + 3/14)0x14, (25) 

q z = ((c S 2 )) r = xi 6 = cr + 3cx 5 r 2 + (-3r + 3ri)cx 7 + (3r - 6r x + 3r 3 )cx 9 

+(-3r 2 + 3r 4 )cxi 3 + (-r + 3r x - 3r 3 + r 5 )cxi 5 , (26) 

where the coefficients in Eqs. (|25|) and (f26| are given as follows 

ro = G(0), 

r a = cosh(JV)G(x)| x= o, 
r 2 = sinh 2 (JV)G(a;)| :I . =0 , 
r 3 = cosh 2 (JV)G(x) 1^=0, 
r 4 = cosh(JV)sinh 2 (JV)G(x)| x=0 , 
r 5 =cosh 3 (JV)G(a;)U = o. (27) 

If we use improved decoupling approximation given in Eq. ([9]) then Eqs. ([23| and (|24|) reduce to the following coupled 
equations 

m = c [cosh( JV) + msinh( JV) + 1 — q z ] 3 F(x)\ x =o, 

= 3cmki + cm 3 fc 3 + (— 6fci + 6k 2 )cmq z + (3fci — 6fc 2 + 3k 4 )cmq 2 z , (28) 



fci = sinh(JV)F(x)| x=0 , 

fc 2 = cosh(JV)sinh(JV).F(x)| x =o, 

fc 3 =sinh 3 (JV)F(a;)| ;c=0 , 

fc 4 = cosh 2 (JV) sinh(JV)F(x)| 2;= o, 



c [gcosh(JV) + msinh(JV) + 1 — q z ] G(x)| x= o, 

cr + 3cm 2 r 2 + (— 3ro + 3r x )cq z + (3r - 6ri + 3r 3 )cg 2 

+ (-3r 2 + 3r 4 )cm 2 q z + (— r + 3n - 3r 3 + r 5 )cg 3 , 



(29) 



where m = ((cjSf)) r and g 
Refs. HHHI [2] ' 



(( Cl (Sf) 2 » r . Eqs. 



and (|2^|) are nothing but just the results obtained in 



44j which exposes the superiority of the present method. 



Now, we should evaluate the thermal and configurational averages of the perimeter site correlations within the 
present formalism. Thus, corresponding to Eqs. (|23|) and (|24[) we have 

(30) 



(C6 ({fs}S z s )) r = (eg ({ft} [co (5«5) 2 cosh( JV) + c S^ sinh( JV) + 1 - c (^) 2 ] ))^ F(x + 7 )| x=0 , 



<c 5 ({/ <5 }(5'|) 2 » r = (c 5 ({/ 5 } [c (^) 2 cosh(JV)+ Co ^sinh(JV) + l- Co (^) 2 ])) r G(x + 7 )U =0 , (31) 



where 7 = (g — represents the effective field produced by the (q — 1) spins outside of the cluster. From Eqs. ([30 
and (|31[) with (5 = {f$} = 1, we can get the perimeter site correlation functions as follows 



((ci;Si)) r = x 4 = aic + a 2 cx\ + (a 3 - ai)cx 



16, 



(32) 



((ciS 2 )) r = x 7 = bic + b 2 cxi + (b 3 - h.)cxie, (33) 

where 

01 = F( 7 ), h = G( 7 ), 

a 2 = sinh( JV)F(x + 7), 6 2 = sinh( JV)G(x + 7), 

a 3 = cosh(JV)F(x + 7 ), fe 3 = cosh( JV)G(x + 7). (34) 

The coefficients in Eqs. (|27|) and (|34|) can be calculated by using the relation e aV /(a;) = /(a; + a). 

For a dilute spin-1 BC model, using Eqs. (1231) . (|2~4"1) . (j3"0)) and (j3"Tj) we derive a set of linear equations by considering 
that (i) the correlations depend only on the distance between the lattice sites, (ii) the average values of a central site 
and its nearest- neighbor site (it is labeled as the perimeter site) are equal to each other with the fact that, in the 
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matrix representations of spin operator S, the spin-1 system has the properties (Sj) 3 = and (S*) 4 = (Sj) 2 - Thus, 
the number of the set of linear equations obtained for the spin-1 Ising system with q — 3 reduces to twenty one, and 
a detailed derivation and the complete set is given inlBl 

Since the thermal and configurational averages of the central site is equal to that of its nearest-neighbor sites within 
the present method then the unknown parameter A in Eq. (|34|) can be numerically determined by the relation 

X\ = X4. (35) 

By solving Eq. (|35|) numerically at a given fixed set of Hamiltonian parameters we obtain the parameter A. Then we 
use the numerical values of A to obtain the site correlation functions such as the longitudinal magnetization ((coSo)) r , 
longitudinal quadrupolar moment ((coS f g)) r and so on, which can be found from Eq. (|B2[) . A = always satisfies 
Eq. (|3"5"1) and gives paramagnetic solution. On the other hand, nonzero solutions of A which satisfy Eq. (|3"5"|) just 
correspond to ferromagnetic state solutions of the system. The critical temperature ksT c / J can be found by solving 
Eq. (133)) in the limit of 7 — > 0. Depending on the Hamiltonian parameters, there may be two solutions [i.e., two critical 
temperature values satisfy Eq. (|35[) ] corresponding to the first (or second) and second-order phase-transition points, 
respectively. We determine the type of the transition by looking at the temperature dependence of magnetization for 
selected values of system parameters. 



III. RESULTS AND DISCUSSION 
A. Site diluted spin-1/2 model 

In Fig. (fTJ) we show the phase diagrams and magnetization, as well as specific heat curves for honeycomb (q = 3) 
and square (q = 4) lattices which can be obtained by solving Eqs. (IA2[) and (|A4[) numerically. In Fig. (JTh) variation 
of magnetization curves are depicted as a function of temperature fc^T '/ J with typical values of site concentration 
c. As expected, we see in Fig. (JTh) that as the temperature increases starting from zero, the magnetization of the 
system decreases continuously, and it falls rapidly to zero at the critical temperature ksTc/J for selected c values. 
The number of interacting sites on the lattice decreases as c decreases and hence, fcgT c / J value of the system and 
the saturation value of magnetization curves also decrease as c decreases. 

In Fig. ([TJd) we examine the effect of site concentration c on the temperature dependence of specific heat of the 
system. We see that as the temperature increases starting from zero, then the specific heat curves exhibit a sharp 
peak at a second-order phase transition temperature which decreases with decreasing c. As c approaches its critical 
value c* at which critical temperature reduces to zero then an additional broad cusp appears and below c* phase 
transition disappears. For c > c* the system forms an infinite cluster of lattice sites however, as c gets closer to 
c* then isolated finite clusters appear and for c < c* the system cannot exhibit long range ferromagnetic order 
even at zero temperature which causes a broad cusp in specific heat vs temperature curves. These observations 
are qualitatively agree with those of Refs. [HI, [ill E EH and show the proper thermodynamic behavior over the 
whole range of temperatures, including the ground-state behavior (C /Nks — >• as ksT/J — > 0) and the thermal 
stability condition (C /Nks > 0). Next, Fig. fT};) represents the variation of the saturation magnetization with site 
concentration. In this figure, we also compare our results (blue line) with those of EFT based on conventional DA 
(C-DA, black line) and improved DA (I-DA, red line) methods. It is clearly evident that site dilution lowers down 
the saturation magnetization. According to C-DA saturation magnetization of the system continuously decreases as 
c decreases then falls rapidly to zero at c* . On the other hand, I-DA predicts a linear decrease at high magnetic atom 
concentrations, but as c decreases gradually then a monotonic decline is observed in the saturation magnetization 
value. On the other hand, according to our results we observe a linear decrement trend up to the vicinity of c* which 
originates as a result of considering the multi-site correlations. Finally, we represent the phase diagram of the system 
in (ksTc/ J — c) plane which separates the ferromagnetic and paramagnetic phases and we compare our results with 
those of the other methods in the literature. According to this figure, critical temperature ksT c / J of system decreases 
gradually, and ferromagnetic region gets narrower as c increases, and A;_bT c / J value depresses to zero at c = c*. Such a 
behavior is an expected fact in dilution problems. Numerical value of critical concentration c* for honeycomb (q = 3) 
and square (q = 4) lattices is given in Table [H and compared with the other works in the literature. It is well known 
that the series expansion (SE) method gives the best approximate values to the known exact results (66|. Therefore, 
we see in Table Q] that the present work improves the results of finite cluster approximation (OSCA and TSCA), as 
well as the other works based on EFT with DA. The reason is due to the fact that, in contrast to the previously 
published works mentioned above, there is no uncontrolled decoupling procedure used for the higher-order correlation 
functions within the present approximation. 
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FIG. 1: Temperature dependence of (a) magnetization, (b) specific heat curves of dilute ferromagnetic system for honeycomb 
(q — 3) and square (q = 4) lattices with some selected values of site concentration c. (c) Ground state magnetizations as a 
function of temperature for q — 3, and q = 4. (d) Phase diagrams of the system in (fesT c /J — c) plane obtained by MFT 
(dash-dotted), DA (dotted), and present work (solid). 



TABLE I: Numerical values of critical site concentration c* for spin-1/2 system obtained within 
the present work for q — 3, 4 and comparison with various approximations in the literature: 
Average coordination number approximation 2/q and Bethe approximation (q — 1) _1 B, RG 
^yjjLCEFT HI, EFT [ll[Il[25|,[Ii, OSCA gaH, TSCA EJEl, CVM [53, MC [33], SE 



q MFT 2/q 


(q-iy 1 RG 


CEFT EFT 


OSCA TSCA CVM MC 


SE 


Present Work 


3 0.667 

4 0.5 


0.5 

0.333 0.602 


0.711 0.5575 
0.558 0.4284 


0.5575 0.5706 0.768 
0.4284 0.4303 0.640 0.413 


0.698 
0.593 


0.6727 
0.4594 



B. Site diluted spin-1 model 



For a site diluted spin-1 BC model defined by Hamiltonian (|20j) . we investigate the thermal and magnetic properties 
of the system by solving Eq. (|B2j) numerically with condition (|35j) . At first, we shall examine the variation of the site 
percolation threshold c* with D/J and ft/ J. In Fig. ([3^) we plot the dependence of the site percolation threshold 
surface with 5.0 < D/J < —1.0 and 0.3<£1/J<2.5. As we can see from Fig. (J2Jl) , the effect of the transverse field 
it/ J on the percolation threshold value clearly depends on the value of the crystal field D/J and vice versa. Namely, 
for the values of ft/ J > 0.565 if we decrease the value of crystal field starting from D/J = 5.0 then c* value increases 
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D/J D/J 

(a) (b) 

FIG. 2: (a) Variation of critical site concentration of a diluted spin-1 BC model for q = 3 with crystal field D/J and transverse 
field Q/J. (b) Critical bond concentration of the same model projected on the same plane. (For interpretation of the references 
to color in these figure legends, the reader is referred to the web version of this article.) 



and and reaches its maximum value. On the other hand, for 0.3 < Q/J < 0.565 c* value increases or decreases 
depending on the value of D/J. Furthermore, for fi/J < 1.56 and sufficiently large positive D/J, c* value remains 
more or less constant and we obtain c* = 0.6727 which is the critical site concentration of spin-1/2 system for q = 3. 
Besides, for D/J = and ft/ J — we get c* = 0.6211 which is higher than the bond percolation threshold value 
of the same system obtained by the same method [Hj]. This value can be compared with the results obtained by 
the other works given in Table [TTJ In Table HH two different critical concentrations obtained by EFT comes from the 
usage of exact or approximate Van der Waerden identity. Using the exact identity one obtains the result of OSCA. By 
comparing Table fl] and Table [TT1 we see that critical site concentration c* of a dilute system depends on the spin value 
S. However, according to the percolation theory (6f| |66| c* only depends on the topology of the lattice and must be 
independent of S. In order to fix this problem, Refs. [27], HH HH suggested to include a positive crystal field D/J 
but, it is clear in Fig. ([2]) that there is an exceptional situation (dark blue region in Fig. (J2Jd)) due to the presence 
of fi/J. Therefore, we can say that topology deformation of the percolation threshold surface illustrated in Fig. @ 
originates from a competition due to the presence of D/J and fi/J in the system. For completeness of the work, we 
also give the critical bond concentration surface of the same model obtained by the same methodology presented in 
this paper for q = 3 (55| . By drawing inspiration from Figs. (J^l) and (03), we think that whether in a site or bond 
dilution problem, the mechanism underlying the complex topological behavior of the critical concentration completely 
originates from a collective effect of both il/J and D/J. 

TABLE II: Site percolation threshold value c* for 
D/J = and fl/J — obtained by present work 
for spin-1 system on a honeycomb lattice. For com- 
parison, the results obtained by OSCA and TSCA 
lH, EFT Hill and SE [gEEl are also S iven - 

OSCA TSCA EFT EFT SE Present Work 
0.5158 0.5449 0.5085 0.5158 0.698 0.6211 



In Fig. ([3]), we represent the phase diagrams of the system in (/esT c /J — D/J) plane for Cl/J = 0, 0.5, 1.0 and 
1.5 where the solid and dotted lines correspond to the second and first order transitions and hollow circles denote 
the tricritical points. The numbers accompanying each curve denote the value of site concentration c. In Fig. ([3]), 
it is obvious that diluting the lattice sites reduces the critical temperatures of the second order phase transitions in 
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FIG. 3: The phase diagrams of a diluted spin-1 BC model in (fc_gT c / J — D/J) plane for q — 3 with selected values of transverse 
field fi/J = 0.0, 0.5, 1.0 and 1.5. Solid and dotted curves correspond to the second and first order phase transitions, respectively. 
Solid circles represent the tricritical points, and the numbers on the curves denote the site concentration c. 



the system for D/J > 0. As seen in the upper left panel of Fig. ([3]), the curve corresponding to pure case (O/J = 
and c = 1.0) exhibits a reentrant behavior of first order where a second order phase transition is followed by a first 
order phase transition at low temperatures for certain negative values of D/J. On the other hand, for c = 0.9 we 
observe an extraordinary feature in the phase diagrams. In other words, there are two regions in D/J plane at 
which a reentrant behavior occurs. The usual one is located within the interval —1.3059 < D/J < —1.0 with a 
tricritical point (D t /J = -1.3058, k B T t /J = 0.4192), and the other is found between -1.0 < D/J < 0.0. The latter 
behavior is quite interesting, since another tricritical point appears at D t /J = —1.0 and ksTt/J = 0.0. Besides, 
for > D/ J > —1.0 the system exhibits a reentrant behavior of second order. On the other hand, if we select c = 0.8 
then we see that the first order phase transitions and tricritical points disappear, and the system exhibits a reentrant 
behavior of second order within the interval —1.1471 < D/J < 0.0. Furthermore, the reentrance disappears as D/J 
becomes positive for all selected values of c. Meanwhile, (ksT c /J — D/J) phase diagrams for some selected values 
of c with n/J — 0.5 are depicted on the upper right panel in Fig. ([3]). It is clearly seen from this figure that the 
system exhibits a first order reentrance for c = 1.0 only in a narrow region —1.4077 < D/J < —1.3796. For c = 0.9, 
tricritical point and reentrance tends to disappear, but if we decrease the magnetic atom concentration further, such 
as for c = 0.8 then the phase diagrams exhibit a bulge with a pronounced second order reentrance within the interval 
—0.9358 < D/J < —0.3554. If we select sufficiently large transverse field strengths, such as fi/J = 1.0, and 1.5 
then the system cannot exhibit first order transitions and tricritical points anymore, even if c = 1.0. In this case, we 
observe only second order phase transitions and ferromagnetic region gets narrower as c decreases. In Ref.[5l|, the 
authors studied the same model for q = 4, but they have not reported the behavior shown in Fig.Q in their paper. 
All of the observations reported here can also be verified by examining the corresponding magnetization curves (see 
FigrU). 

The evolution of the phase diagrams shown in Fig. ^ for fl/J — and 0.5 are depicted in Fig. (Q}. As seen 
in Fig. (QJi) where il/J = 0, the phase diagrams exhibit a reentrant behavior of second order within the interval 
— 1.0 < D/J < 0.0 for c = 0.80, 0.78, 0.76 and 0.75 whereas the reentrant phase transition region for c = 0.74 and 0.73 
is divided into two parts: The first part is located in the vicinity of D/J = —1.0 which gets narrower as c decreases 
while the second one is observed between —0.3607 < D / J < 0.0. If we decrease c further, such as for c = 0.69 and 0.68 
(see Fig. (HJd)) reentrance disappears. However, for 0.62 < c < 0.67 another reentrant regime appears, but now for 
D/J > which gets narrower as c decreases. Similarly, Figs. ((4}:) and (J4ji) represents the evolution of phase diagrams 
of the system corresponding to the upper right panel in Fig. © where fi/J = 0.5. As seen in Fig. (0J:), the system 
undergoes only a second order phase transition for c = 0.87. However for lower site concentrations e.g. c = 0.83 and 
0.81, reentrant behavior of second order appears between —0.609 < D/J < —0.4938 and —0.7891 < D/J < —0.386, 



11 



0.9 




-1.2 -0.9 -0.6 -0.3 0.0 0.3 0.6 0.9 0.0 0.5 1.0 1.5 2.0 2.5 



D/J D/J 

(c) (d) 

FIG. 4: Evolution of the phase diagrams given in Fig. © for (a) fi/J = 0.0 and 0.73 < c < 0.80, (b) Q./J = 0.0 and 
0.62 < c < 0.69, (c) Q/J = 0.5 and 0.72 < c < 0.87, (d) Q/J = 0.5 and 0.63 < c < 0.69. The numbers accompanying each 
curve denote the site concentration c. 

respectively. In addition, for c < 0.79 the system exhibits a bulge in the reentrant phase transition regime which gets 
narrower as c decreases. For c < 0.72 (Fig.((4}l)), reentrance disappears, and the system undergoes a second order 
phase transition form a paramagnetic state to a ferromagnetic state with increasing temperature. For c < 0.63, site 
concentration of the system approaches to percolation threshold c*, hence the ferromagnetic region becomes fairly 
narrow. 

Next, as a complementary investigation of Fig. (J3]), we investigate the effects of transverse field interactions ft/ J 
on the phase diagrams of the system in Fig. (0 for some selected values of c. The upper left panel in Fig. §5§ 
corresponds to the phase diagrams of the pure system [67| where we see that reentrant phase transitions tend to 
disappear and tricritical points decrease as Jl/J decreases for c = 1.0. Moreover, according to the upper right panel in 
Fig. ([5]), another ferromagnetic phase boundary arises between —1.0 < D/J < 0.0 and —0.833 < D/J < —0.0434 for 
c = 0.9 and fi/J = 0.0, 0.25, respectively. These additional phase transition lines get narrower as O/J increases, and 
disappear after a certain value of il/J. Furthermore, phase diagrams shown in lower left and right panels invariably 
exhibit second order phase transitions for c = 0.8 and 0.7 which is independent from transverse field value, and the 
reentrance is not observed anymore for c = 0.7 and Q/ J > 0.5. 

In Fig. (J6]), we examine the phase diagrams of the system in a (fcsT c /J — c) plane for D/J = —1.25, —0.5, 0.0 and 
1.0 with some typical £l/J values. As seen in this figure, the system exhibits a tricritical behavior and its critical 
temperature cannot reach zero for D/J = —1.25 and Q/J = 0.0, hence we cannot speak on any percolation threshold 
value. However, phase transition temperature of the system reduces to zero at c* = 0.9764 for il/J — 0.25. It is 
clear from upper left panel in Fig. © that c* value for il/ J = 0.75 is greater than those of SI/ J = 0.5, however 
it is lower than those of fl/J = 0.25. On the other hand, for D/J = —1.25 and Q./J — 0.5 the first order phase 
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FIG. 5: The phase diagrams of a diluted spin-1 BC model in (fesTc/J — D/J) plane for q = 3 with selected values of site 
concentration c = 1.0,0.9,0.8 and 0.7. Solid and dotted curves correspond to the second and first order phase transitions, 
respectively. Solid circles represent the tricritical points, and the numbers on the curves denote the value of transverse field 
Q/J. 



transitions and tricritical points disappear, and we barely observe a second order reentrance which also disappears 
for fl/J = 0.75. If we select D/J = —0.5 then we cannot see any evidence of first order phase transitions, and critical 
site concentration c* of the system decreases as il/J increases up to £1/ J = 0.75. For il/J> 0.75 we observe that 
c* value tends to increase as fl/J increases which is consistent with the indications depicted in Fig. (|2Jl) . Similar 
discussions are also valid for D/J = 0.0 and 1.0. Additionally, as an interesting characteristic of the system, we may 
note that the conditions for the occurrence of a second order reentrance in the system is rather complicated, since the 
existence or extinction of reentrance is rather sensitive to the collective effects of D/J, Vi/ J and c. 

As a final investigation, let us represent the temperature dependence of the magnetization for some selected values 
of Hamiltonian parameters corresponding to the phase diagrams depicted throughout Figs. In Fig. ([7]), 

the typical transition profiles are shown for several values of Hamiltonian parameters. For example, if we select 
D/J = —1.25 and fl/J = 0.0 magnetization curves corresponding to c — 1.0, 0.95 and 0.90 exhibit a discontinuous 
jump at a first order transition temperature then gradually decrease and reduce to zero at a second order phase 
transition temperature with increasing temperature. This is an example of reentrance of first order in which a first 
order transition is followed by a second order transition. As an example of second order reentrant behavior, we 
can take a look at the magnetization curves that exist in a second order reentrant regime. For instance, we see 
that the magnetization curves exhibit two critical temperatures of the second order for c = 0.90, 0.85 and 0.80 with 
D / J — —0.5 and Q/ J = 0.0. On the other hand, as an example of a second order ferromagnetic-paramagnetic phase 
transition, magnetization curves exhibit two different characteristics. As an example of the first case, we see that 
as the temperature increases then the magnetization of the system falls gradually from its saturation magnetization 
value at fc^T/J = 0.0 and decreases continuously up to the vicinity of the transition temperature and vanishes at a 
critical temperature ksT c / J for c = 1.0 with D/J = —0.5 and il/J — 0.0 (corresponding to pure case), whereas in the 
second case the magnetization of the system exhibits a temperature-induced maximum with increasing temperature 
which is depicted on the lower left and right panels of Fig. (J7J . 



IV. CONCLUDING REMARKS 



In this work, we have investigated the thermal and magnetic properties of a site diluted spin- 1/2 Ising model 
and a spin-1 Blume Capel (BC) model in the presence of transverse field interactions. We have introduced an 
effective-field approximation that takes into account the multi-site correlations in the cluster of a considered lattice 
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FIG. 6: The phase diagrams of a diluted spin-1 BC model in (fcsT c /J — c) plane for q — 3 with selected values of the crystal 
field D/J — —1.25,-0.5,0.0 and 0.5. Solid and dotted curves correspond to the second and first order phase transitions, 
respectively. Solid circle denotes the tricritical point, and each curve in panels is plotted for a specific transverse field value. 
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FIG. 7: Temperature dependence of the magnetization curves as functions of Hamiltonian parameters D/J, Q/J and c. Solid 
and dotted magnetization curves exhibit second and first order phase transition properties, respectively. 



with an improved configurational averaging technique. Our method is capable of locating the possible first order 
phase transition temperatures, as well as tricritical points, and under certain simplifications, equations of state 
obtained within the present approximation can be reduced to those obtained by conventional or improved decoupling 
approximation techniques which exposes the superiority of the present work. 

For a spin- 1/2 Ising system, we have obtained results that are superior to those estimated by conventional mean 
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field theory (MFT) and effective field theory (EFT) based on a decoupling approximation, especially for the critical 
site concentration (i.e. site percolation threshold) value c* for honeycomb (q = 3) and square (q — 3) lattices. Our 
estimated values c* = 0.6727 and c* = 0.4594 for q = 3 and q = 4, respectively are the best approximate values to 
the results of MC and SE methods among the other works based on MFT or EFT. 

In particular, we have investigated the phase diagrams and magnetization curves of a site diluted spin-1 BC model 
in the presence of transverse field interactions and we have shown that diluting the lattice sites may cause some drastic 
changes on some of the characteristic features of the model. For this model, we have examined the variation of the site 
percolation threshold c* with the crystal and transverse field interactions which has not been reported in the literature 
before. In the absence of crystal and transverse fields, the percolation threshold value of a site diluted spin-1 model 
for q = 3 is estimated as c*=0.6211 which improves the results obtained by other EFT based approximations. In 
addition, we have found that the percolation threshold value c* strictly depends on the value of crystal and transverse 
field interactions, as well as the topology of the lattice. We have also given the global phase diagrams, especially the 
first order phase transition lines that include reentrant phase transition regions. The results presented in this paper 
clearly indicate that the conditions for the occurrence of a second order reentrance in the system is rather complicated, 
since the existence or extinction of reentrance is rather sensitive to the competing effects between D/J, fl/J and c. 
These observations cannot be observed by ignoring any of these Hamiltonian parameters in the system. 

As a result, we can conclude that all of the points mentioned above show that our method improves the conventional 
EFT methods based on decoupling approximation. Therefore, we hope that the results obtained in this work may be 
beneficial from both theoretical and experimental points of view. 
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Appendix A: Derivation of complete set of linear equations for spin-1/2 Ising Model 

In the present formalism, all of the site correlations including central, as well as perimeter site magnetizations are 
denoted by For instance, for a honeycomb lattice (q = 3) we have 

Xi = ((c 5 )) r , 

x 2 = ((coSqCiSi)},., 

%3 = ((coSoCiSiC 2 S2))r, 

x 4 = ((ciSi)) r , 

X5 = ((ciSlC 2 S 2 }}r, 

X6 = ((ciSlC 2 S 2 C 3 S 3 )} r - (Al) 

Basis correlation functions for central and perimeter sites are defined respectively as follows: 

m = ((coSo)) r = x± = (3c — 6c 2 + 3c 3 )x 4 Ki + (6c 2 — 6c 3 )x 4 K 2 + Zc^XiK?, + cxqK 4 , 
{{ciSi)) r = Xi = {c-c 2 )A 1 +c 2 A 2 + cx 1 A 3 . 

By expanding Eq. ([B]) with i — and {fi} = c 4 Si we get 

((ciS , ic S'o))r = x 2 = (3c 2 - 6c 3 + 3c 4 )Ki + (6c 3 - 6c 4 )K 2 + 3c 4 K 3 + c 2 x 5 K 4 . 

By the same way, putting i — and {/,;} = ciS\c 2 S 2 in Eq. (|6]) we obtain 

x 3 = (-3c 2 + 3c 3 ):E4^i + (6c 2 - 6c 3 )x 4 K 2 + 3c 3 x 4 K 3 + c 3 x 4 K 4 . 

Similarly, by using Eq. (fTB")) with 6=1 and {fs} = c 2 S 2 , and {fs} = c 2 S 2 c 3 S 3 we find x$ and Xq, respectively. Hence, 
we get the complete set of correlation functions as follows: 



X\ 


= (3c -6c 2 4 


- 3c 3 )x 4 K 4 + (6c 2 - 6c 3 )x 4 K 2 + 3c 3 x 4 K 3 + cx 6 K> 


x 2 


= (3c 2 -6c 3 


■f 3c 4 )K 4 + (6c 3 - 6c 4 )K 2 + 3c 4 K 3 + c 2 x 5 K 4l 


x 3 


= (-3c 2 + 3c 


3 )x 4 K 1 + (6c 2 - 6c 3 )x 4 K 2 + 3c 3 x 4 K 3 + c 3 x 4 K 4l 


x 4 


= (c-c 2 )^ 


+ c 2 A 2 + cxxA 3 , 


x 5 


= cx 2 A 3 , 




xe 


= cx 3 A 3 . 





(A2) 



where the coefficients K n , (n = 1, 4) and Aj, (I = 1, 3) are given in Eqs. (|13[) and (fTB)) . respectively. 
On the other hand, corresponding to Eq. (|A1|) . for a square lattice (q = 4) we have 



Xi = ((c S )) r , 

X 2 = ((c Q S CiSi)) r , 

x 3 = {{c S ciSic 2 S 2 }} r , 

x 4 = {{c S ciSic 2 S 2 c 3 S 3 }} r , 

x 5 = ((ciSi)) r , 

xe = {{ciSic 2 S 2 )} r , 

x 7 = ((ciSic 2 S 2 c 3 S 3 )) ri 

x s = ((ciSic 2 S 2 c 3 S 3 c 4 S 4 )) r . (A3) 

By following the same procedure given for q = 3 above, we get the complete set of linear equations for a square 
lattice (q = 4) as follows: 

x 4 = (4c - 12c 2 + 12c 3 - 4c 4 )x 5 L 1 + (12c 2 - 24c 3 + 12c 4 )x 5 L 2 + (12c 3 - 12c 4 )x 5 L 3 

+Ac 4 x 5 L 4 + (4c - 4c 2 )x 7 L 5 + 4c 2 x 7 L 6 , 
x 2 = (4c 2 - 12c 3 + 12c 4 - 4c 5 )L! + (12c 3 - 24c 4 + 12c 5 )i 2 + (12c 4 - 12c 5 )L 3 + 4c 5 L 4 

+ (4c 2 - 4c 3 )x 6 L b + 4c 3 x 6 L 6 , 
x 3 = (-8c 2 + 12c 3 - 4c 4 )x 5 L 1 + (12c 2 - 24c 3 + 12c 4 ):E5.L2 + (12c 3 - 12c 4 )x 5 L 3 
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where 



X4 

x 5 
x e 
x 7 



+4c 4 X5L 4 + (4c 3 - 4c i )x 5 L 5 + A^x^Lq, 
{Ac 2 - 4c 3 )a; 6 ii + (-12c 2 + 12c 3 )x 6 L 2 + (12c 2 
+(4c 2 - 4c 3 )a; 6 L5 + Ac 3 x 6 L 6 , 
(c - c 2 )Bi + c 2 B 2 + CX1B3, 
(c - c 2 )x 5 Bi + c 2 x 5 B 2 + cx 2 B 3 , 



12c 3 )x 6 L 3 + Ac 3 x 6 L 4 



c 2 x 5 B 2 
(c — c 2 )x 6 B 1 + c 2 x e B 2 
c 2 x 7 B 2 



{c-c 2 )x 7 B l 



cx 3 B 3 , 
cx 4 B 3 , 



L 2 

L 3 - 

L 6 



■■ sinh(JV) taiih.(/9a;)| 3 , = o, 

: cosh(JV) sinh( JV) tanh(/3a;)| x= o, 

cosh 2 ( JV) sinh(JV) tanh(/3a:)| x= o, 

cosh 3 ( JV) sinh( JV) tanh(^x)\ x=0 , 

sinh 3 ( JV) ta,nh((3x)\ x=0 , 

cosh(JV) sinh 3 (JV) tooh{Px)\ x=0 . 



B 2 
B 3 



tanh(^(a; + 7))| x=0 , 
cosh(JV) tanh(/3(x + 7))| x=0 , 
sinh( JV) tanh(/3(x + ^))\ x=0 , 



(A4) 



Phase diagrams and magnetization curves can be obtained by solving Eq. 

Xi = x 5 . 



4]) numerically with the condition 



(A5) 



Appendix B: Derivation of complete set of linear equations for spin-1 BC model 

We label the site correlations as i — 1, 2, 21. The complete list is as follows 





X\ = ({coS )) r , 


X12 


= ((c S ciSlc 2 Sl))n 




x 2 = ((c Q SQCiSi)) r: 


Xl3 


= ((ciSiC 2 S 2 C 3 S 3 }}r, 


X3 


= {{coSociSic 2 S 2 )) r , 


Xl3 


= ((ciS , ic 2 S , |c 3 S'3)) I ., 




X 4 = ((ciSi)) r , 


X15 


= ((cxSfaSlcsS 2 )),., 




%5 = ((ciSiC 2 S 2 )) r , 


Xl6 






= ((ciSlC 2 S 2 C 3 S 3 )} r , 


Xyj 


= ((coS%CiSl}) r , 




X 7 = ((cxSf^r, 


X\% 


= {{c SQCiS 2 }} r , 




X 8 = {{ciS\C 2 Sl)) r , 


Xig 


= ((c SQClSlC 2 S 2 ))r, 




XQ = {{c\S 2 C 2 S 2 )) r , 


X 2 Q 


= ((c SQClSlC 2 S 2 )) r , 




xw = ((coS ciS 2 )) r , 


X21 


= {{c SlciSlc 2 Sl)) ri 


Xll 


= ((coS ciSic 2 S 2 )) r . 







The correlation functions Xi, i = 1,2,3 are obtained from Eq. (|23|) . For example, putting i = and {/;} = c\S\ 
and {fi} — ciS\c 2 S 2 in Eq. (|23[) we obtain x 2 and x 3 correlation functions, respectively as follows 

((coiSoCiS'i))r — x 2 — 3ck\x 7 + (— 6fci + 6k 2 )cx g + ck 3 xi 3 + (3ki - 6k 2 + 3k 4 )cxi 5 , 
{{coSoCiSic 2 S 2 )) r = x 3 — (— 3fci + 6fc 2 )cx 8 + (3/ci - 6fc 2 + fc 3 + 3fc 4 )cxi4, 

The equations labeled Xj with j — 4,5,6 are derived from Eq. (1301) . In a similar way, the correlation functions Xk 
with k = 7, 8, 15 and xi with I — 16, 17, 21 can be easily obtained by using Eqs. ([24| and (f3l"|) . respectively. By 
following the above procedure, we can get the complete set of linear equations as follows: 



x\ — 3cx 4 k\ + cx§k 3 + (— 6fci + 6fc 2 )cxs + (3fci — 6fc 2 + 3^4)0^14, 

x 2 = 3ck\x 7 4- (— 6fci + 6fc 2 )cxg + ck 3 x\ 3 + (3fcx — 6fc 2 + 3fc4)ca;i5, 

x 3 = (— 3fei + 6k 2 )cxs + (3ki — 6k 2 + k 3 + 3fc4)cxi4, 

x 4 = die + a 2 cxi + (a 3 - ai)cx 16 , 



17 



x 5 = a x cx 4 + a 2 cx 2 + (a 3 - a\)cxn, 

x 6 = (L1CX5 + a 2 cx 3 + (a 3 - ai)cxi9, 

X7 = &ic + fo 2 cxi + (63 - ^i)cxi 6 , 

x$ = b 1 cx 4 + b 2 cx 2 + (63 - h)cxi 7 , 

Xg = biCX 7 + b 2 cx 10 + (63 - &i)cxi 8 , 
xio = b 2 cx 16 + b 3 cxi, 
xu = b 2 cx 17 + b 3 cx 2 , 
X12 = b 2 cxis + b 3 cx w , 
X13 = bicx 5 + b 2 cx 3 + (63 - bi)cxw, 
xu = bicx s + b 2 cxn + (63 - bi)cx 20 , 
xi 5 = bxcx 9 + b 2 cxi 2 + (63 - 61)0x21, 

xie = cr + 3cr 2 x 5 + (-3r + 3ri)cx 7 + (3r - 6ri + 3r 3 )x 9 

+(-3r 2 + 3r 4 )cx 13 + (-r + 3r x - 3r 3 + r 5 )cx 15 , 
x 17 = (-2r + 3rx)cx 4 + (3r + 3r 2 + 3r 3 - 6rx)cx 8 

+(-r + 3ri - 3r 2 - 3r 3 + 3r 4 + r 5 )cxi 4 , 
xis = (-2r + 3ri)cx 7 + (3r - 6ri + 3r 2 + 3r 3 )cx 9 

+(-r + 3ri - 3r 2 - 3r 3 + 3r 4 + r 5 )xi 5 , 
X19 = (r - 3ri + 3r 2 + 3r 3 )cx 5 

+(-ro + 3n - 3r 2 - 3r 3 + 3r 4 + r 5 )cxi 3 , 
x 2 o = (ro - 3n + 3r 2 + 3r 3 )cx 8 

+(-r + 3ri - 3r 2 - 3r 3 + 3r 4 + r 5 )cxi 4 , 
x 2 i = (r - 3ri + 3r 2 + 3r 3 )cx 9 

+(-r + 3ri - 3r 2 - 3r 3 + 3r 4 + r 5 )cxi 5 . (B2) 
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